#################
# Figure 2: Long term change in y after a change in x at t=1-
#################

rm(list=ls())
y<-c(0,1)
x<-c(1)
i<-2

while(round(y[i], 5)!=round(y[i-1],5)){
  i<-i+1
  y[i]<-0.2*y[i-1]+0.5*x+0.3*x*y[i-1]+0.5*1
  print(y[i], quote=F)

}




rm(list=ls())
y<-c(1.99,1.99,1.99,1.99,1.99,1.99,1.99,1.99,1.99,2)
x<-c(1.5)
i<-10

while(round(y[i], 5)!=round(y[i-1],5)){
  i<-i+1
  y[i]<-0.2*y[i-1]+0.5*x+0.3*x*y[i-1]+0.5*1
  print(y[i], quote=F)

}

t<-1:i
plot(y~t, type="l", xlab="time")
abline(v=10, lty=2)

dev.copy2eps(file="evolutionofY.eps")




#################
# Figure 1: The (Instantaneous) Effect of x_it on y_it
#################

rm(list=ls())
require(mvtnorm)
set.seed(123809)
x<-runif(30*10, min=-2, max=2)
z<-runif(30*10, min=-2, max=2)
y<-c()
unit.ind<-c()                                       
time.ind<-c()
lag.y<-c()

# generating the data set 

for(j in 1:30){
  for(i in 1:10){
  
    time.ind[10*(j-1)+i]<-i
    unit.ind[10*(j-1)+i]<-j
    
    if(i!=1){
      y[10*(j-1)+i]<-0.2*y[10*(j-1)+i-1]+0.5*x[10*(j-1)+i]+0.3*x[10*(j-1)+i]*y[10*(j-1)+i-1]+0.5*z[10*(j-1)+i]+rnorm(1,mean=0,sd=0.4)
      lag.y[10*(j-1)+i]<-y[10*(j-1)+i-1]
    }else{
      y[10*(j-1)+i]<-runif(1, min=-1, max=1)
      lag.y[10*(j-1)+i]<-NA
    }
  
  }
}


# running the model 

lagTx<-lag.y*x
model<-lm(y~lag.y+lagTx+x+z)
summary(model)

# drawing the coefficients out of the asymptotic VCV
beta.draws<-rmvnorm(1000, mean=coefficients(model), sigma=vcov(model))

# collect and plot the MEs for each value of y_lag
y.seq<-seq(from=-6, to=4, by=0.1)
me<-matrix(data=NA, ncol=3, nrow=length(y.seq))

for(i in 1:length(y.seq)){

  y.test<-y.seq[i]
  me.draws<-beta.draws[,4]+beta.draws[,3]*y.test
  me[i,]<-quantile(me.draws, probs=c(0.025, 0.5, 0.975))
  

}

plot(me[,2]~y.seq, type="l", ylim=c(-2, 2), xlim=c(-6,4), main="Instantaneous Marginal Effect", xlab=expression(y[it-1]), ylab=expression(partialdiff*y[it]/partialdiff*x[it]))
lines(me[,1]~y.seq, lty=2)
lines(me[,3]~y.seq, lty=2)


dev.copy2eps(file="instantME.eps")




#################
# Figure 3: The Long Term Effect of x_it on y, from x=0 to x=x_new
#################

# calculate long term marginal effects and 95% CIs out of the beta draws from before
x.new<-seq(from=-1.5, to=1.5, by=0.05)
ss.hi<-c()
ss.lo<-c()
ss.med<-c()
for(i in 1:length(x.new)){

  x.lo<-0
  x.hi<-x.new[i]
  ss.draws<-((beta.draws[,4]*x.hi+beta.draws[,5]*0.5)/(1-beta.draws[,2]-beta.draws[,3]*x.hi))-((beta.draws[,4]*x.lo+beta.draws[,5]*0.5)/(1-beta.draws[,2]-beta.draws[,3]*x.lo))
  ss.quant<-quantile(ss.draws, probs=c(0.025, 0.5, 0.975))
  ss.hi[i]<-ss.quant[3]
  ss.lo[i]<-ss.quant[1]
  ss.med[i]<-ss.quant[2]

}


plot(ss.med~x.new, type="l", ylim=c(-1,4.2), main=expression(paste("Long Term Marginal Effect of Change from ", x==0," to ", x[new])), xlab=expression(x[new]), ylab=expression(paste("(", y, "|", z==0.5,", ", x==x[new], ") - (",y, "|", z==0.5,", ", x==0, ")")))
lines(ss.hi~x.new, lty=2)
lines(ss.lo~x.new, lty=2)
legend("bottomright", lty=c(1,2), legend=c("marginal effect", "95% CI"))

dev.copy2eps(file="longtermME.eps")



#################
# Figure 4: The Long Term Effect of z_it on y, plotted against x_it
#################

# calculate long term marginal effects and 95% CIs out of the beta draws from before
x.new<-seq(from=-1.5, to=1.5, by=0.05)
ss.hi<-c()
ss.lo<-c()
ss.med<-c()
for(i in 1:length(x.new)){

  x.lev<-x.new[i]
  ss.draws<-((beta.draws[,4]*x.lev+beta.draws[,5]*1)/(1-beta.draws[,2]-beta.draws[,3]*x.lev))-((beta.draws[,4]*x.lev+beta.draws[,5]*0)/(1-beta.draws[,2]-beta.draws[,3]*x.lev))
  ss.quant<-quantile(ss.draws, probs=c(0.025, 0.5, 0.975))
  ss.hi[i]<-ss.quant[3]
  ss.lo[i]<-ss.quant[1]
  ss.med[i]<-ss.quant[2]

}


plot(ss.med~x.new, type="l", ylim=c(0.25,2.55), main=expression(paste("Long Term Marginal Effect of ", Delta[z]==1)), xlab=expression(x[o]), ylab=expression(paste("(", y, "|", z==1,", ", x==x[o], ") - (",y, "|", z==0,", ", x==x[o], ")")))
lines(ss.hi~x.new, lty=2)
lines(ss.lo~x.new, lty=2)
legend("bottomright", lty=c(1,2), legend=c("marginal effect", "95% CI"))

dev.copy2eps(file="longtermME_z.eps")

